#ifndef _BLIP3_H_
#define _BLIP3_H_

//#include "Equations/TwoPhaseEquation.hpp"
#include "./TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class BLIP3;

  template <typename T>
  struct Typeinfo <BLIP3<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <BLIP3<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <BLIP3<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

  
  template<typename Tprec, int Dim>
  class BLIP3 : public TwoPhaseEquation<BLIP3<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<BLIP3<Tprec, Dim> > TP_BLIP3;
    
    using GeneralEquation< TP_BLIP3 >::aE;
    using GeneralEquation< TP_BLIP3 >::aW;
    using GeneralEquation< TP_BLIP3 >::aN;
    using GeneralEquation< TP_BLIP3 >::aS;
    using GeneralEquation< TP_BLIP3 >::aB;
    using GeneralEquation< TP_BLIP3 >::aP;
    using GeneralEquation< TP_BLIP3 >::sp;
    using GeneralEquation< TP_BLIP3 >::dx;
    using GeneralEquation< TP_BLIP3 >::dy;
    using GeneralEquation< TP_BLIP3 >::dz;
    using GeneralEquation< TP_BLIP3 >::bi;
    using GeneralEquation< TP_BLIP3 >::ei;
    using GeneralEquation< TP_BLIP3 >::bj;
    using GeneralEquation< TP_BLIP3 >::ej;
    using GeneralEquation< TP_BLIP3 >::bk;
    using GeneralEquation< TP_BLIP3 >::ek;
    using GeneralEquation< TP_BLIP3 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_BLIP3 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_BLIP3 >::applyBoundaryConditions3D;

    using TP_BLIP3::S;
    using TP_BLIP3::phi_0;
    using TP_BLIP3::Srw;
    using TP_BLIP3::Sro;
    using TP_BLIP3::mu_w;
    using TP_BLIP3::mu_o;
    using TP_BLIP3::k;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    BLIP3() : TwoPhaseEquation<BLIP3<prec_t, Dim > >() { }    
    ~BLIP3() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D() { };
    inline bool calcCoefficients3D() { };
    inline void printInfo() { std::cout << " BLIP3 "; }
  };

/*
 *  Lineal for realtive permeability , Average for Sw
 */
template<typename Tprec, int Dim>
inline bool BLIP3<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;

    // Lineal 
    static prec_t mult_o = k / (1 - Srw - Sro) * mu_o * dx ;
    static prec_t mult_w = k / (1 - Srw - Sro) * mu_w * dx ;
    
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {      

      // Average
      Sw_e = ( S(i+1) + S(i) ) * 0.5;
      Sw_w = ( S(i-1) + S(i) ) * 0.5;      

      // Lineal
      aE (i) = (1 - Sro - Sw_e) * mult_o + (Sw_e - Srw) * mult_w ;
      aW (i) = (1 - Sro - Sw_w) * mult_o + (Sw_w - Srw) * mult_w ;
      aP (i) = aE (i) + aW (i);      
    }
    applyBoundaryConditions1D();
    return 0;
}



} // Tuna namespace


#endif //_BLIP3_H_

















